Social network dynamics of face-to-face interactions 
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The recent availability of data describing social networks is changing our understanding of the 
"microscopic structure" of a social tie. A social tie indeed is an aggregated outcome of many 
social interactions such as face-to-face conversations or phone calls. Analysis of data on face-to- 
face interactions shows that such events, as many other human activities, are bursty, with very 
heterogeneous durations. In this paper we present a model for social interactions at short time 
scales, aimed at describing contexts such as conference venues in which individuals interact in small 
groups. We present a detailed analytical and numerical study of the model's dynamical properties, 
and show that it reproduces important features of empirical data. The model allows for many 
generalizations toward an increasingly realistic description of social interactions. In particular in 
this paper we investigate the case where the agents have intrinsic heterogeneities in their social 
behavior, or where dynamic variations of the local number of individuals are included. Finally we 
propose this model as a very flexible framework to investigate how dynamical processes unfold in 
social networks. 
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I. INTRODUCTION 

In the last decade complexity theory has greatly ad- 
vanced thanks to the availability of extensive data on 
a wide variety of networked systems. Empirical stud- 
ies have uncovered the presence of ubiquitous features 
in complex networks, such as the small- world property, 
or strong heterogeneities in the topological structure, re- 
vealed for instance by broad degree distributions [lJ-Q] ■ 
These findings have deeply affected our understanding of 
self-organized networks and have been used for the in- 
vestigation, characterization and modeling of many dif- 
ferent systems such as infrastructure or biological and 
social networks. 

Many works have studied the influence of complex net- 
work topologies observed in real networks on the dynam- 
ical phenomena that unfold on them [(| 0] ■ While these 
studies have mostly focused on networks considered as 
static objects with a fixed topology, networks' structures 
may in principle evolve, links may appear and disappear. 
A first approach consists of assuming that links are cre- 
ated and annihilated at a constant rate, independently of 
the dynamical process that takes place on them |8l— flOj] - 
Networks can however display more interesting properties 
such as an adaptative behavior, in which the dynamics 
on the network and of the network are related by feed- 
back effects (TT|-[r7j. Moreover, empirical investigations 
have shown that link durations can significantly deviate 
from a Poisson process [l8l - [23j . 

Social networks 0, [H[ are prominent examples of 
evolving networks. Social relationships are indeed con- 
tinuously changing, possibly in a way correlated with the 
dynamical processes taking place during social interac- 
tions (such as epidemic spreading or opinion dynamics). 
Consequently, a number of works have been devoted to 



modeling the dynamics of social interactions. Issues in- 
vestigated in this context are in particular community 
formation [26l - l28j and the evolution of adaptive dynam- 
ics of opinions and social ties through schematic models 
in which links can disappear or be rewired at random 
[TTMT7I . Moreover, social networks evolve on many dif- 
ferent timescales. The static representation of social ties 
hide indeed dynamical sequences of events such as face- 
to-face interactions, phone calls or email exchanges, and 
can be measured by aggregating fast social interactions 
over a certain period of time. 

Recently, technological advances have made possible 
the access to data sets that give new insights into such 
link internal dynamics, characterized by sequences of 
events of different durations. Traces of human behav- 
ior are often unwittingly recorded in a variety of con- 
texts (financial transactions, phone calls, mobility pat- 
terns, purchases using credit cards, etc.). Data have 
been gathered and analyzed about the mobility patterns 
inside a city [29| , between cities f30j, a s well as at the 
country and at worldwide levels [20, I3ll - l33j . At a more 
detailed level, mobile devices such as cell phones make it 
possible to investigate individual mobility patterns and 
their predictability [HI, [35[ . Mobile devices and wearable 
sensors using Bluetooth and Wifi technologies give ac- 
cess to proximity patterns of pairs of individuals jl8l l36l — 
[39||, and even face-to- face presence can be resolved with 
high spatial and temporal resolution [2l| - [23l l40j . Fi- 
nally, on-line interactions occurring between individuals 
can be monitored by logging instant messaging or email 
exchange [4lT - [47l ]. 

The combination of these technological advances and of 
heterogeneous data sources allows researchers to gather 
longitudinal data that have been traditionally scarce in 
social network analysis (48|,|4!| ■ Analysis of such data sets 
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has clearly shown the bursty nature of many human and 
social activities, revealing the inadequacy of many tra- 
ditional frameworks that posit Poisson distributed pro- 
cesses. In particular, the durations of "contacts" between 
individuals, as defined by the proximity of these indi- 
viduals, display broad distributions, as well as the time 
intervals between successive contacts [HI, [III I21M23I l39j . 
Burstiness of interactions has strong consequences on dy- 
namical processes [H, l50l - l54j |. and should therefore be 
correctly taken into account when modeling the inter- 
action networks. New frameworks are therefore needed, 
which integrate the bursty character of human interac- 
tions and behaviors into dynamic network models. 

While a lot of modeling efforts have been devoted to 
static networks, the development of models of dynamic 
networks has indeed until recently attracted less atten- 
tion P3, tn m m m. in & recent p & p er & we 

have presented an agent-based modeling framework to 
describe how individuals interact at short times scales in 
venues such as social gatherings (e.g., scientific confer- 
ences). The model is based on a reinforcement dynamics 
(in the spirit of the preferential attachment in complex 
networks 57J and of Hebbian learning) which might be 
responsible for the bursty dynamics of social face-to-face 
interactions. The proposed mechanism implies that the 
longer an agent is interacting in a group, the smaller 
is the probability that he/she will leave the group; the 
longer an agent is isolated the smaller is the probability 
that he/she will form a new group. 

In the present paper, we present an extensive char- 
acterization of the model's properties, and show that it 
reproduces important features of empirical data on so- 
cial interactions at short timescales. We characterize the 
rich phase diagram of the model, which includes station- 
ary and non-stationary regions. The analysis of the dy- 
namical properties of the model shows that it yields sta- 
tionary broad distribution of group lifetimes even if the 
underlying dynamics is non-stationary, as also found in 
empirical data. In order to illustrate the model's versa- 
tility, we give two examples of how it can be extended 
to more realistic cases: in the first example, agents can 
have heterogeneous propensities to form group with oth- 
ers; in the second example, we introduce the possibility 
of a varying population, where the number of individuals 
can be an arbitrary function, or extracted from empiri- 
cal data. In the two cases we show how properties very 
close to the ones of real- world data sets can be obtained. 
The proposed model is easily implementable, uses simple 
but realistic mechanisms, reproduces a certain number of 
empirical facts, and is amenable to further refinements. 
It can therefore be used to produce artificial data sets 
of bursty interaction networks on which dynamical phe- 
nomena can be simulated and studied. 

The paper is organized as follows. In Sec. HH we re- 
view the main properties of a representative empirical 
data set describing the face-to-face proximity of individ- 
uals in social gatherings. Section. Mil is devoted to the 
definition of the modeling framework and to the analyt- 



ical and numerical study of the simplest versions of the 
model. Section. IIVI is devoted to two extensions of the 
model. We outline some conclusions and perspectives in 
Sec. H 



II. EMPIRICAL DATA 

The infrastructure developed by the SocioPatterns 
project, described in (2TI !2il. has yielded measurements 
about the face-to-face proximity of individuals in differ- 
ent types of social contexts (hospital, primary school, 
scientific conferences, museum), with a fine grained time 
resolution. The infrastructure is currently based on radio 
frequency identification devices (RFID). Individuals par- 
ticipating to the data collection are asked to wear small 
RFID tags on their chests (as a conference badge), that 
emit radio packets at very low power. The parameters of 
the infrastructure are tuned so that the tags can exchange 
radio packets only when the individuals wearing them 
face each other at close range (about 1 to 1.5 m), and 
so that face-to-face proximity events ("contacts") can be 
assessed with very high accuracy with a time resolution 
of 20s. We present here for the sake of completeness 
some characteristics of the data collected during the 6th 
European Semantic Web Conference (ESWC, Heraklion, 
Greece) in 2009. Analysis of other data sets can be found 
in Refs. (2l| - [23j . with very similar features. 

Face-to-face proximity patterns have been collected for 
175 voluntary participants (among the 305 conference at- 
tendants) over 3 days. 14520 contact events have been 
registered, corresponding to 27.3 contacts per individual 
per day. The contact durations display a very broad dis- 
tribution, close to a power law, (shown in Fig. [TJ: the 
average duration is of 46 s, but long-lasting contacts are 
as well observed, and no characteristic contact timescale 
can be extracted from the data. Figure [T] also displays 
a quantity of interest, namely the distribution of time 
intervals between the start of two consecutive contacts 
of a given individual A with two distinct persons B and 
C . In other words, if A starts a contact with B at time 
tAB ' and then starts a different contact with C at £ac> 
the inter-contact interval is defined as r = tAc — tAB- 
These time intervals constrain causal processes such as 
information diffusion or epidemic spreading, as they de- 
termine the timescale after which an individual receiving 
some information or disease is able to propagate it to an- 
other individual. As shown in Fig. [TJ broad distributions 
are also found in this case. We also show in Fig. [2] the 
distributions of lifetimes of groups of size n + 1 (n = 
corresponds to an isolated person, n = 1 to a pair, etc). 
All these distributions are broad, compatible with power- 
law shapes, and become narrower for increasing n (larger 
groups are less stable than smaller ones). 

The burstiness of the contact pattern revealed by the 
broad distribution of contact durations has consequences 
also on aggregated views of the dynamical contact net- 
work. Aggregated contact networks over a given time 
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FIG. 1: Distribution of contact durations (a) and of time 
intervals between two contacts involving a common individual 
(b). 
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FIG. 2: (Color online) Distributions P n (j) of the durations 
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window are denned as follows: each node corresponds to 
an individual, and an edge is drawn between two nodes 
if at least one contact event has been registered during 
the time window between the two corresponding individ- 
uals. Each edge is weighted by the sum of the contact 
durations between these individuals during this time win- 
dow. As shown in Fig. [3] (top right), the distributions of 
such weights are broad, independent of the time window 
considered (see also [22, HH ) • 

Figure [3] also displays other characteristics of aggre- 
gated networks constructed with time windows of dif- 
ferent lengths. As also found in [22, HH for other de- 
ployments of the same infrastructure, the distribution of 
degrees (the degree gives the number of distinct individ- 
uals with whom a given individual has been in contact) 
are not broadly distributed. This behavior is in con- 
trast with the degree distribution of many empirically 
studied social networks @, HH, [59[ . It should however be 
emphasized that we are here dealing with face-to-face in- 
teractions occurring in a restricted environment among a 



FIG. 3: (Color online) Aggregated networks' characteristics 
for the empirical ESWC data: (a) degree distribution, (b) 
weight distribution, (c) average strength of nodes of degree k, 
vs k, and (d) average Herfindahl-Hirschman index of nodes of 
degree k, vs k. 

relatively small pop ulation and on short timescales, while 
studies such as |58l . l59j are concerned with social ties cre- 
ated and defined over much longer times cales. Moreover, 
we note that a narrow degree distribution (a power law 
with a large exponent) has also been found in the degree 
distribution of networks defined by phone-calls data sets 

The strength Si of a node i is defined as the sum of the 
weights Wij between i and its neighbors j, that is, 

s i =^Tw ij . (1) 

3 

It gives the cumulated durations of the interactions of the 
corresponding individual. The average strength of nodes 
of degree k [(s(k))] indicates how the weights are dis- 
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tributcd. If the weights are uniformly distributed among 
the links of the networks, the average strength (s(k)) 
grows linearly with k, that is, (s(k)) oc k(w). On the 
contrary, if stronger ties are more frequently linked to 
highly connected agents, a superlinear behavior of (s(fc)) 
versus k is observed JUJ . The RFID data on face-to face 
interaction are consistent with a linear or slightly faster 
behavior of (s(k)) versus k [22|, [23|, hinting at a weak cor- 
relation between weights and degrees. Correlations be- 
tween network topology and distribution of the weights 
have been found in various complex networks with broad 
degree distribution [3l|, and several models based on re- 
inforcement dynamics have been proposed in this context 

Another measure of the weights' distribution is given 
by the Herfindahl-Hirschman index Y% [63], [64|, also 
known in the physics literature under the name of "par- 
ticipation ratio" . This index, defined as 

where Af(i) refers to the set of neighbors of i, gives a 
measure of the heterogeneity of the weights among the 
neighbors of a node. When all weights Wij of the links 
connected to node i are equal, that is, Wij = Si/ki, this 
index is inversely proportional to the degree fcj of node i, 
that is, 

Y 2 {i) = I (3) 

On the contrary, when one of the links has a much larger 
weight than the others, Y 2 (i) is close to 1. The departure 
of kiY 2 (i) from 1 thus indicates the local heterogeneity 
of weights around each node. Figure [3] shows that the 
average of this quantity over nodes of degree k [kY 2 (k)j, 
is larger than 1 and increases with k: each individual di- 
vides his/her time unevenly among her/his contacts, and 
this behavior becomes more pronounced as the number 
of distinct contacted persons increases. 

III. A MODEL FOR SOCIAL DYNAMICS AT 
SHORT TIMESCALES 

In this section, we present a model describing the social 
dynamics of N agents forming small groups of different 
size. In this model, each agent can interact with any 
other agent. This model describes therefore social inter- 
actions within a closed environment of relatively small 
size where agents are free to meet, such as a conference 
venue. We assign to each agent i = 1, 2, . . . , N a coordi- 
nation number rii = 0, 1, 2 . . . , N indicating the number 
of agents interacting with him/her. If an agent i has 
coordination number rii = he/she is isolated, and if 
rii = n > he/she is part of a group of n + 1 agents, 
who all interact with each other (thus forming a clique). 
We also assign to each agent i the temporal variable U 



indicating the last time at which his/her coordination 
number rii has changed. 

The dynamics of the model is as follows. Starting from 
random initial conditions, at each time step t the follow- 
ing steps are performed: 

(1) An agent i is chosen randomly. 

(2) The agent i updates his/her coordination number 
rij = n with a certain probability p n (t, U) that may 
depend on the agent's state, on the present time t, 
and on the last time tj at which i's state evolved. 
With probability 1 — p n (t, U), the agent does not 
change state. 

If the coordination number m is updated, the ac- 
tion of the agent is chosen with the following rules. 

(i) If the agent i is isolated, that is, rii = 0, 
he/she starts an interaction with another iso- 
lated agent j chosen with probability propor- 
tional to pa(t, tj). The coordination number of 
the agent i and of the agent j are then updated 
according to the rule rii — > 1 and rij — > 1. 

(ii) If the agent i is interacting in a group, that 
is, rii — n > 0, with probability A the agent 
leaves the group and with probability (1 — A) 
he/she introduces an isolated agent in the 
group. If the agent i leaves the group, his/her 
coordination number is updated (n, — > 0), as 
well as the coordination numbers of all the 
agents in the original group, that is, — > 
n — 1 (for all agents k in the original group). 
On the contrary, if the agent i introduces an- 
other isolated agent j to the group, j is chosen 
with probability proportional to po(t,tj) and 
the coordination numbers of all the interact- 
ing agents are changed according to the rules 
rij — > n + 1 and — > n + 1 (for all k in the 
group). 

The structure and properties of the interactions be- 
tween agents depend on the choice of the probabilities 
p n , which control the tendency of the agents to change 
their state, and on the parameter A, which determines 
the tendency either to leave groups or on the contrary to 
make them grow. The simplest choice consists of consid- 
ering constant probabilities p n (t,t') = p n : at each time, 
every agent has a fixed probability to form a group or 
split from a group. In this case the formation of the 
groups is a Poisson process, and the distributions of the 
durations of contacts between agents, or the lifetime of a 
group, are exponentially distributed. 

As recalled in the introduction however, the distribu- 
tions of the times describin g h uman activities are typi- 
cally broad 0, [H, HI M, SI, HI HI, and are clearly 
closer to power-laws that lack a characteristic time scale 
than to exponentials. A possible explanation of such re- 
sults is given by mechanisms in which the decisions of 
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the agents to form or leave a group are driven by mem- 
ory effects dictated by reinforcement dynamics, that can 
be summarized in the following statement: the longer an 
agent is interacting in a group, the smaller is the proba- 
bility that he/she will leave the group; the longer an agent 
is isolated, the smaller is the probability that he/she will 
form a new group. In particular, such reinforcement prin- 
ciple implies that the probabilities p n (t, t') that an agent 
with coordination number n changes his/her state de- 
pend on the time elapsed since his/her last change of 
state. A simple way to introduce this hypothesis is to 
consider functions p n (t,t') = p n (t — t'). Reinforcement 
mechanisms are then described by decreasing functions 
p n . We will see in the next subsections how the evolu- 
tion equations of the number of agents in each state can 
be written at the mean-field level for arbitrary functions 
p n . Finding solutions of this set of equations is however 
not always possible, and we will focus on functions p n 
scaling as l/(t— t 1 ), for two reasons: on the one hand, 
it represents one of the cases that is fully amenable to 
analytical computations and, on the other hand, such a 
scaling behavior is needed to obtain power-law distribu- 
tions for the contact durations and thus dynamical char- 
acteristics compatible with empirical data. We consider 
in particular functions given by 



Pn(t,t') 



l + (t-t')/N' 



(4) 



and, in order to reduce the number of parameters, we 
moreover take b n = b\ for every n > 1, indicating the fact 
that interacting agents change their state independently 
on the number of other agents n with whom they are 
interacting, provided that n > 1. The model's parameter 
are thus bo, b\, and A. 



A. Pairwise interactions 

Let us first consider a restricted version of the model, 
in which the agents can only interact in pairs. This set- 
up is obtained by setting A = 1 and by considering initial 
conditions in which the agents interact at most in groups 
of size 2. In this case, each agent is thus assigned a 
variable ni = 0, 1 indicating if the agent i is isolated 
{iii = 0) or interacting with another agent (n^ = 1). 

As in the analysis of empirical data, the most imme- 
diate quantities of interest concern the time spent by 
agents in each state, the duration of contacts between 
two agents, and the time intervals between successive 
contacts of an agent. To gain insight into these temporal 
properties of the system, we can write rate equations for 
the evolution of the numbers N n (t,t') of agents in state 
n = at time t who have not changed state since time 
t' . In the mean-field approximation, and treating time 
and numbers as continuous variables, these equations are 



given by 

dN (t,t') Mt,t') 



dt 



N 



P0 (t,t')+TT 1Q (t)5tt> 



= -2 — — Pl{t,t ) + 7T01 (*)<)«/, (5) 



dt 



N 



where the transition rates 7r„ im (i) denote the average 
number of agents switching their states from n to m 
(n — > m) at time t. If the agents make their decisions 
according to the reinforcement dynamics described by 
the probabilities p n (t,t') given by Eq. (j4|), the dynamic 
equations ([5]) have a solution of the form 



N (t,t') = w w (t') 1 + 



JVi(f,f) = 7r i(f) 1 + 



t-t' 
N 

t-t' 
N 



-2b 



-2&i 



(6) 



Since the total number of isolated agents who change 
their state at time t is equal to ttqi (t) and the total num- 
ber of interacting agents who change their state is equal 
to 7Tio(t), it follows that 7Tio(i) and 7Toi(i) are given in 
terms of N (t,t') and Ni(t,t') by the relations 



r=i 
2 * 

n °i® = ^ &o(M')^o(M')- 
t'=i 



(7) 



To solve the coupled set of equations ^ and (J7J, we 
assume self-consistently that ttiq (t) and ttqi (t) are either 
constant or decaying in time as power laws. Therefore, 
we assume 



7i"io(*) = ttio 



7I"0l(*) = ttqi 



N 



N 



(8) 



To check the self-consistent assumption Eq. (JSJ, we 
insert it in Eqs. ([6]) and and compute the values of 
the parameters ao, Oix, tti and 7r i that determine the 
solution in the asymptotic limit t — > oo. If ao = ai = 0, 
we obtain a stationary solution in which 7Tio(t) = 7Tio and 
ftoiit) = ^oi i arc independent of time. On the contrary 
if ao > or ai > 0, the system is non-stationary, with 
transition rates 7Tio(t) and 7Toi(t) decaying in time. The 
system dynamics slows down. In Appendix \K\ we give 
the details of this self-consistent calculation in the large 
N limit, which yields ao = a.\ — a and tt\o = 7Toi = tTj 
with 

a = max(0, 1 — 26i, 1 — 26 ) , 
sin [2ir min (60, bl 



[1 -<*(«,<>)] 



(2b - l)(26x - 1) 
2(6o + 6i-l) 



S(a,0). 



(9) 
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FIG. 4: (Color online) Phase diagram of the pairwise model. 
The white area indicates the stationary regime in which the 
transition rate is constant. The colored (gray) area indicates 
the non-stationary phase. 
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FIG. 5: (Color online) Evolution of the transition rate 7Tio(t) 
in the different phase regions of the pairwise interaction 
model. The simulation is performed with N = 1000 agents 
for a number of time steps T max = N x 10 , and averaged over 
10 realizations. The simulations are performed in the station- 
ary region with parameter values bo = bi =0.7 (circles) and 
in the non-stationary region with parameter values bo — 0.3, 
b\ = 0.7 (squares) and bo = bi =0.1 (triangles). The lines 
indicate the analytical predictions Eqs. 



FIG. 6: (Color online) Probability distribution of the du- 
rations of contacts Pi(t) and of the inter-contact durations 
Pq(t) in the stationary region, for the pairwise model. The 
data is reported for a simulation with N = 1000 agents, run 
for T max — N x 10 5 elementary time steps, with parameter 
values bo = 0.6, bi = 0.8. The data is averaged over 10 real- 
izations. 



Eq. ((9j) . In this stationary state the number of iso- 
lated agents and the number or interacting agents 
are constant on average, but the dynamics is not 
frozen, since w > 0: agents continuously form and 
leave pairs. The simulations shown in Fig. [5] for 
bo = b\ = 0.7 confirm this analytical prediction. 



• Non stationary region (bo < 0.5 or b\ < 0.5J - In 
this region of the phase diagram, the self-consistent 
equation predicts a non-stationary solution with 
7Tio(t) and 7Toi(f) decaying with t as a power-law 
of exponent a = max(l — 2bo,l — 26i). Figure [5] 
shows such a decay for bo = 0.3, b\ = 0.7 and for 
bo =b\ =0.1, which is however truncated by finitc- 
size effects for t larger than t c (N) oc N. There- 
fore the system eventually becomes stationary with 
a very slow dynamics (very small transition rates 
7Tio(£) and 7T i(i)). 



The analytically predicted dynamical behavior or the 
model can be summarized by the phase diagram depicted 
in Fig. 0] (that we discuss now in more detail), together 
with the numerical simulations of the stochastic model 
displayed in Fig. O 

• Stationary region (bo > 0.5 and b\ > 0.5j - In 
this region of the phase diagram, the self-consistent 
equation predicts a = 0, so that a stationary state 
solution is expected, where 7Toi(i) = ft is given by 



Empirical studies often focus on the statistics of con- 
tact durations between individuals, and of the time in- 
tervals between two contacts of a given individual. These 
quantities of interest can be computed in our model, re- 
spectively, as the probabilities -Pl(t) that an agent re- 
mains in a pair during a time r = (t — t')/N, and 
Pq(t) that an agent remains isolated for a time interval 
t = (t — t')/N. These probabilities are determined by the 
numbers of agents in each state and the rates at which 
the agents change their state. The probability distribu- 
tions of the time spent in each state, integrated between 
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FIG. 7: (Color online) Probability distribution of the du- 
rations of contacts -Pi(r) and of the inter-contact durations 
Po (r) in the non- stationary region of the pairwise model, with 
bo < 0.5 and bi < 0.5. In this region we observe some devia- 
tions of the probabilities Pi(t) and Po(t) from the power-law 
behavior for large durations. The data are reported for a 
simulation with N = 1000 agents run for T max = N x 10° 
elementary time steps, with parameter values bo = 6i = 0.1. 
The data are averaged over 10 realizations. 



the initial time and an arbitrary time t, are given by 

rt-Nr 



P„(r) cx 



p n (t' + Nr,f)N n (f + N T ,t')dt' (10) 



for n = 0, 1. Inserting the expression given by Eq. ^ 
for N n (t,f) and the definition of p n (t,t') given by Eq. 
(H|) in Eq. (JTUJ), we obtain the power- law distributions 



P„(r) <x(1+t) 



-2h„-l 



(11) 



for n = 0, 1. These analytical predictions are compared 
with numerical simulations in Fig. [5] for bo = 0.6, b\ = 
0.8 (stationary system) and in Fig. [7] for b = b\ = 0.1 
(non-stationary ttiq and 7Toi). Interestingly, even when 
the system is non-stationary, the distributions P n (r) re- 
main stationary. 



B. Formation of groups of any size 

In this subsection we extend the solution obtained for 
the pairwise model to the general model with arbitrary 
value of the parameter A, where groups of any size can 
be formed. Therefore the coordination number rii of each 
agent i can take any value up to N — 1. Extending the 
formalism used in the previous subsection, we denote by 
N n (t, f) the number of agents with coordination number 
n = 0, 1, ...,N — 1 at time t, who have not changed 
state since time f . In the mean field approximation, the 
evolution equations for N n {t,f) are given by 



dN (t,f) 
dt 

dt 

dN n (t,f) 
dt 



N (t,t') 



N 
' JV 



P1 {t,t') 



i>l 



+[no A (t)+Tr 2 ,i(t)}5tt>, 

N n (t,t') , „ 
-(n+1) "y ' piitj) 



+[ir n - lti (t)+ir n+1 , n (t)+Tr 0>n (t)}5tt>, n>2. (12) 

In these equations, the parameter e(t) indicates the 
rate at which isolated nodes are introduced by another 
agent in already existing groups of interacting agents. 
Moreover, n mn (t) indicates the transition rate at which 
agents change coordination number from m to n (i.e. 
m — > n) at time t. In the mean-field approximation the 
value of e(t) can be expressed in terms of N n (t,t') as 



e(t) 



E n >iELi^(M>i(MO 

El=iNo(t,t')p (t,t>) 



(13) 



In the case of reinforcement dynamics described by the 
probabilities p n {t,t') given by Eq. Q, and assuming 
that asymptotically in time e(t) converges to a time- 
independent variable, that is, lim t _ i . 0O e(t) = e, the so- 
lution to the rate equations f| T2|) in the large time limit 
is given by 



N (t,t') 

N n (t,t') 
with 



N (t',t') 
N n (t',t'] 



t-f 
N 

t-f 
N 

t-f 
N 



-6 [2+(l-A)e] 



-26, 



-(n+l)6i 



(14) 
for n > 2, 



N (f,t') = Y,n n ,o(f), 



n>l 



N n (f,f) 



7T ,l(O +^2,1(0, (15) 

(t r ) +TT0,n(t') for n > 2. 



The transition rates n m ,n(t) can be determined in 
terms of N n (t,f) as shown in the Appendix IB"1 In order 
to solve the equations we make the further assumption 
that the transition rates n mn (t) are either constant or 
decaying with time according to a power law, that is. 



N 



(16) 



Self-consistent calculations (see Appendix [Bj) determine 
the value of the quantities e, a mn , and TT mn . For A > 0.5 
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FIG. 8: (Color online) Phase diagram of the general model 
with formation of groups of arbitrary size. The region behind 
the green surface corresponds to the stationary phase [i.e., 
region (I), with A > 0.5, 61 > 0.5 and b > §jE±]. The region 
in front of the green surface and above the blue one [region 
(II)] corresponds to a non-stationary system with decaying 
transition rates. Strong finite size effects with a temporary 
formation of a large cluster are observed in the region below 
the blue surface [i.e., region (III) with A < 0.5]. 



the self-consistent assumption Eq. <|16[) is valid and we 
find, as in the case of pairwise interactions, that a m , n = a 
V(m, n), with 



3A — 1 

a = max ( 0, 1 - b Q ^ - ^ , 1 - 2&i 



(17) 



This solution generalizes the case of the pairwise model, 
which is recovered by setting A = 1. For A < 0.5 the self- 
consistent assumption breaks down and we will resort to 
numerical simulations. 

The probability distributions of the time spent in each 
state, integrated between the initial time and an arbi- 
trary time t, are given by 



Pn(r) cc 



t-Ni 



p n (t' + NT,t')N n (t' + NT,t')dt' . (18) 



Inserting the expression given by Eq. ([Til) for N n (t,t') 
and the definition of p n (t,t') given by Eq. ((4J in Eq. 
(|18p , we obtain the power-law distributions 



P (t) oc (1+r) 



-6o[2+(l-Ae)]-l 



P„(t) oc (l + r)"^ 1 ) 61 - 1 forn>l. 



(19) 



The phase diagram of the model is summarized in Fig. 
|U We can distinguish between three phases. 

• Region (I) -The stationary region: b\ > 0.5, 60 > 
(2A - 1)/(3A - 1) and A > 0.5- In this region, the self- 
consistent solution yields a — 0: the transition rates 



FIG. 9: (Color online) Transition rate 7rio(t) for the model 
in the presence of groups of any size, for different parameters 
A, 60, 61 corresponding to the different regions of the phase 
diagram. The straight lines correspond to the analytical pre- 
dictions. The simulation is performed with N = 1000 agents 
for a number of time steps T max = N X 10 4 . The data are 
averaged over 10 realizations. 
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FIG. 10: (Color online) Distribution -Pn(r) of durations of 
groups of size n + 1 in the stationary region. The simulation 
is performed with TV = 1000 agents for a number of time steps 
Tmax — N x 10 5 . The parameter used are bo — b\ — 0.7, A = 
0.8. The data are averaged over 10 realizations. The dashed 
lines correspond to the analytical predictions Eqs. (|19|) . 



7Tmn(i) converge rapidly to a constant value (see Fig. 
|9] for a comparison between numerics and analytics for 
7Tio(£)) and the system reaches a stationary state. In 
Fig.[TU]we compare the analytical solution given by Eqs. 
(|19p with the numerical simulations in the stability re- 
gion, finding perfect agreement. As predicted by Eqs. 
(|T9l . P n {j) decays faster as n increases: larger groups are 
less stable than smaller ones, as found in the empirical 
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FIG. 11: (Color online) Distribution of time intervals between 
successive contacts of an individual, for A = 0.8, bo = 0.7 and 
bi = 0.3 and 0.9. The simulation is performed with N = 10 4 
for a number of time steps T max = N x 10°. The data are 
averaged over 10 realizations. 



FIG. 13: (Color online) Distribution Pn(r) of durations of 
groups of size n + 1 in the non stationary region [region (II)] . 
The simulation is performed with N = 1000 agents for a num- 
ber of time steps Tmax = N x 10 5 . The parameter used are 
bo — 0.3 and b\ — 0.7, A = 0.8. The data are averaged over 
10 realizations. The dashed lines correspond to the analytical 
predictions Eqs. (|19[l . 




FIG. 12: (Color online) Average coordination number (n) vs 
A for bo — bi — 0.7 . The simulation is performed with 
N = 2000 agents for a number of time steps T max = N x 10 3 . 
(n) is computed in the final state over 30 realizations. The 
solid line indicates the theoretical prediction given by Eq. 
(2D). 



data sets. Figure QT] displays the distribution Pab-ac( t ) 
of time intervals between the start of two consecutive con- 
tacts of a given individual, which is as well stationary and 
displays a power-law behavior. 

The average coordination number (n) is given by 



(n) 



7T10 n ( n + 1) 



2A ^ (n+l)bi - 1 



1 - A 
A 



(20) 



where the detailed calculation and the value of nio(t) 
are given in Appendix [B] This expression diverges as 
A ^ 0.5. In Fig. [T^] we show the perfect agreement 
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FIG. 14: (Color online) Average coordination number (n) as 
a function of time in the Region II of the phase diagram for 
different values of the parameters A, bo and b\ . The data is in 
very good agreement with the theoretical expectations given 
by Eqs. (|21|l — (22). The simulations are performed with 
N = 1000 agents for a number of time steps T maa; = N x 10 4 . 
The data are averaged over 10 realizations. 



between the result of numerical simulations of (n) and 
the theoretical prediction. 

• Region (II) -Non- stationary region: b\ < 0.5 or 
b < (2A - 1)/(3A - 1), and A > 0.5 - The dynamics 
in this region is non-stationary and the transition rate 
is decaying with time as a power-law, as shown in Fig. 
[5] where we report ir w (t) as a function of t. Neverthe- 
less, the distributions of lifetimes of groups of various 
sizes P n (r), and of inter-contact times Pab-ac{ t ), rc_ 
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FIG. 15: (Color online) Average coordination number (n) 
for A = 0.2, bo = &i = 0.7. The simulations of a single 
realization are performed with N = 250 and N — 500 agents, 
respectively, for a number of time steps T max = N x 10 . 



main stationary. These distributions are shown in Figs. 
1 131 and [Til In this region, the average coordination num- 
ber in the limit t/N ^> 1 remains small, even as A — > 0.5. 
In particular from the mean-field solution of the dynam- 
ics (see appendix [Bj) the theoretical solution of the model 
predicts that, for A > 0.5 and t — > oo 



and 



(n) = 1 for a 



(n) = for a = 1 



2&i, 



3A- 1 
1 2A - 1 ■ 



(21) 



(22) 



Figure [Ml shows the agreement of this predicted behavior 
with simulation results for several values of bo an d bi and 
A = 0.7. In this region, as A — > 0.5 with fixed bo and bi, 
we have a = 1 — 2bi and (n) — > 1. Therefore no diverging 
behavior is observed. 

• Region (III) Strong dependence on the number of 
agents N and non- stationary dynamics: A < 0.5 - In 
this region the self-consistent assumption given by Eq. 
(|16[) breaks down, and we find numerically that the av- 
erage coordination number (n) strongly depends on the 
number of agents N and on time. In order to give a typ- 
ical example of the corresponding dynamical behavior, 
Fig. [15] displays (n) as a function of time for two single 
realizations of the model corresponding to two different 
values of N. Interestingly, the distributions of lifetimes 
of groups of various sizes P n {i~) remain stationary even 
in this parameter region (not shown). 



C. Aggregated networks 

In the previous paragraphs, we have shown how our 
modeling framework produces dynamical properties of 



the interactions between agents that yield broad distri- 
butions of contact and inter-contact times. In order to 
understand the structure of the resulting interaction net- 
works at coarser temporal resolutions, it is as well in- 
teresting to investigate the properties of the aggregated 
networks, constructed as in Sec. [TTJ Given a starting 
time to and a temporal window AT, the nodes of these 
networks are the agents, and a link is drawn between 
two agents whenever they have been in contact between 
to and to + AT, with a link weight given by the total 
time during which they have interacted in [to , to + AT] . 
As in Sec. QTJ the degree fcj of an agent i is given by 
the number of distinct agents with whom i has been in 
contact in [to, to + AT], while its strength Sj is the sum 
of the interaction times with other agents, and the par- 
ticipation ratio Yi (i) quantifies the heterogeneity of the 
times spent by i with these other agents. 

As an exhaustive exploration of the aggregated net- 
works and of how their properties depend on the model's 
parameter would be tedious, we simply report in Fig. QJj] 
the properties of aggregated networks for increasing win- 
dow lengths AT and for two sets of parameters. Some 
properties are qualitatively similar to the empirically ob- 
served networks. In particular, the degree distributions 
are peaked around an average value that increases with 
AT. As time passes, each agent encounters more and 
more distinct other agents, and the distribution P(k) 
globally shifts towards larger degrees. The links weights 
distributions are broad, and extend to larger values as 
AT increases. Some other properties seem to depend 
strongly on the model's parameters. In particular, the 
average strength of nodes of degree k, and the average 
participation ratio of nodes of degree k, can have shapes 
rather different from the empirical ones. Moreover, the 
time window lengths AT on which the aggregated net- 
work remains sparse are rather restricted. 



IV. EXTENSIONS OF THE MODEL 

A. Heterogeneous model 

In the previous section, we have assumed that all the 
agents have the same tendency to form a group or to leave 
a group, that is, the probabilities p n do not depend on the 
agent who performs a status update. Real social systems 
display however additional complexity since the social 
behavior of individuals may vary significantly across the 
population. A natural extension of the model presented 
above consists therefore of making the probabilities p n 
dependent on the agent who is updating his/her state. 
To this aim, we assign to each agent i a parameter r\i 
that characterizes his/her propensity to form social in- 
teractions. In real networks this propensity will depend 
on the features of the agents [65||. In the model we as- 
sume that this propensity, that we call "sociability", is 
a quenched random variable, which is assigned to each 
agent at the start of the dynamical evolution and re- 
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who have been interacting since time t'. The mean-field 
equations for the model are then given by 

dN (t,t',r,) n N (t,t', V ) 

— m — = 2 n "ob*'") 

+ 7r io(*)< 5 tf , 
dAW, 77,77') Nx{t,t' \r,,r,') 



at 



N 



[pi(*.* / >'?)+Pi(<>*W)] 
(24) 
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With the expression for p n (t,t',r]) given by Eqs. (|23| we 
find 

t-t'\ -2i 



N (t,f,r,) = nUt')(l+- 1 f-) 



10 w io _ 



10 



(h) 



10 



/ t — t'\ -2+?7+»?' 

N^tJ^rf) = + — ) • (25) 

The transition rate 7r^ gives the rate at which agents 

with r/i G [77, 77 + A77] become isolated, and ir^ is the rate 
at which pairs ij with r/i G [rj, 77 + A77], rjj G [77', 77' + A77] 
are formed. These rates can be expressed as a function 
of No(t, t', 77) and N±(t, t' , 77, 77') according to 
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FIG. 16: (Color online) Aggregated networks' characteristics , 
for the model with constant number of agents (N = 250), ^oi CO = 2 / 
for time windows AT of increasing lengths, and two sets t'.t" 
of parameters: (60, 61, A) = (0.55,0.8,0.9) (a, b, c, d) and 
(0.7,0.7,0.8) (e, f, g, h). 



AW, V,rf) 



N 



Pi (*>*'>'/)]> (26) 



C{t)N 



Po(t,t',v)Po(t,t",v'), 



mains constant, and we assume for simplicity that it is 
uniformly distributed in [0,1]. In this modified model, 
the probability p l n (t,t') that an agent i with coordina- 
tion number n since time t' changes his/her coordination 
number at time t is given by 



Po(M') = 
K(M') = 



where C(t) is a normalization factor given by 
t 

C{ 1) = E E v)po(t, t', r,). (27) 

t' = l Tj 

To solve this problem with the same strategy used for the 
model without heterogeneity we make the self-consistent 
assumption that the transition rates are either constant 
or decaying as a power law with time: 



l + {t-t')/N' 
l + (t-t')/N' 



for n > 1. 



(23) 



In this setup, the parameters (bo,bi), which did not de- 
pend on i in Eq. are replaced for each agent i by 
the values (77^, 1 — 77^): a large 77^ corresponds to an agent 
who prefers not to be isolated. 

The agents' heterogeneity adds a significant amount of 
complexity to the problem, and we have reached an ana- 
lytical solution of the evolution equations only in the case 
of pairwise interactions (A = 1). The general case can be 
studied through numerical simulations as we discuss at 
the end of this section. 

Let us denote by No(t,t',rj) the number of isolated 
agents with parameter 77^ G [77, 77 + A77] who have not 
changed their state since time t' . Similarly, we indicate 
by N±(t, t' , 77, 77') the number of agents in a pair joining 
two agents i and j with 77, G [77, 77 + A77] , rjj G [?/ , 7/ + A77] , 



t \ - a (v,v') 



(28) 
(29) 



In appendix [C] we give the details of the self-consistent 
calculation, which leads to the analytical prediction 



(30) 



a{rj) = max I 1 - 2rj, 77 - - j , 
a(i7,r/') = a(i]) + a(rf) , 



and the value of 7r^ is given by 



~ri _ ) B(l-2r),2n) 



7T 



r?<i 



10 



B(v-h,l) 11 ~ 2 



(31) 



In order to check the validity of our mean-field calcu- 
lation, we study the probability distribution Pq(t) of the 
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FIG. 17: (Color online) Distributions of times spent in state 
and 1 for the heterogeneous model. The simulation is per- 
formed with iV = 10 4 for a number of time steps T max = 
N x 10 5 . The data are averaged over 10 realizations. The 
symbols represent the simulation results (circles for n = 
and squares for n = 1). The dashed lines represent our ana- 
lytical prediction. In order to improve the readability of the 
figure we have multiplied Pi(r) by a factor of 10 _1 . 



FIG. 18: (Color online) Distribution P^{r) of contact dura- 
tions of individuals with sociability r\ in the pairwise heteroge- 
neous model. The simulations are performed with N = 1000 
agents and Tmax = N x 10 5 time steps. The data are aver- 
aged over 10 realizations. The data decays as a power-law 
Pi(t) oc t~^( v \ and we report the exponents £(rf) as a func- 
tion of r\ in the inset. 



durations of inter-contact periods and the distribution 
P\{t) of the durations of pairwise contacts, which are 
given, when averaged for a total simulation time T max , 
by 

pT maI -NT si 

P (r) oc / dt ^(^(l+r)- 2 "- 1 , 



Pi(t) cc 



Jo 

T max -Nr f l r l 



dt I dt] I dr/nW 



(2-t ? - ? /)(1 + t)'"+"'- ; \ 



(32) 



where p(?y) is the probability distribution of 77. In Fig. [T7] 
we compare the probabilities of intcrcontact time Pq(t) 
and contact time Pi (r) averaged over the full population 
together with the numerical solution of the stochastic 
model, showing a perfect agreement. In Fig. [15] more- 
over, we show the distributions Pi{t) of the contact du- 
rations of agents with rji £ (77, 77 + A?7). Power- law be- 
haviors are obtained even at fixed sociability, and the 
broadness of the contact duration distribution of an agent 
increases with the" sociability" of the agent under consid- 
eration. 

As previously mentioned, the model can be extended 
by allowing the formation of large groups, by setting 
A < 1. The results of numerical simulations performed 
for a particular value of A are shown in Fig |19l Power 
law distributions of the lifetime of groups are again found 
and, as in the basic model without heterogeneity of the 
agents, larger groups are more unstable than smaller 
groups, as Pn(j~) decays faster if the coordination num- 
ber n is larger. As the parameter A — > 0.5 there is a 
phase transition and the average coordination number 
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FIG. 19: (Color online) Distribution P n (j) of the durations of 
groups of size n+1 in the heterogeneous model with formation 
of groups of any size. The data are shown for simulations of 
iV = 1000 agents performed = N x 10 5 time steps 

and A = 0.8, averaged over 10 realizations. 



diverges. In Fig. [20] we show that (n) oc (A — 0.5) -<5 
with S = 1 within the statistical fitting error, similarly 
to what happens in the homogeneous case. Overall, the 
main features of the model are therefore robust with re- 
spect to the introduction of heterogeneity in the agents' 
individual behavior. 
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sidcrcd 

• the agents who leave the system can re-enter it: 
the global number of agents (absent and present) 
is constant. We will see that a set-up in which an 
agent who has left the premises cannot re-enter the 
system leads to an interesting modification of the 
system's properties. 
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FIG. 20: (Color online) Average coordination number (n) of 
the agents in the heterogeneous case as a function of A. The 
solid line indicates the best fit with (n) oc (A - 0.5) ~ s with 
S = 0.996 in agreement with the exponent —1 within the 
statistical uncertainty. The data correspond to simulations of 
TV = 500 agents performed over T m ax = TV x 10 3 time steps. 
The data are averaged over 10 realizations. 



B. Model with variable number of agents 

When measurements on the proximity or face-to- 
face contacts of individuals are performed, the number 
of agents present on the premises typically fluctuates. 
Moreover, the activity of the agents fluctuates as well, 
because for instance of day /night patterns or periods of 
coffee/ lunch breaks in a conference. This is in particular 
the case in the data sets described in Sec. [TTJ 

In this subsection, we use these simple remarks to put 
forward an extension of our model of interacting agents 
that mimics in a more realistic way the data gathering 
process, and that can be used to produce more realistic 
artificial data sets. The main point is to consider, instead 
of a population with a fixed number of agents TV, an open 
system leaving the possibility for agents to leave or enter 
it. To this aim, we simply introduce for each agent the 
possibility to be in a state called "absent" . In this "ab- 
sent" state, agents may be isolated or not, but the mea- 
surement infrastructure is not able to know it. Statistics 
on the time spent by an agent in a given state, or of the 
duration of contact and inter-contact times, are thus ob- 
tained by considering only "present" agents. Overall, the 
model's dynamics is as described in Sec. IIII1 with param- 
eters bo, &i, A, with the addition that at each time step, 
"absent" agents can enter the system, or agents can leave 
the system and become "absent" . Various rules could be 
thought of, but for the sake of simplicity we consider that 
an agent i entering the system become isolated (n, ; — > 0), 
and that agents of any state n can leave the system (one 
could also allow "absent" agents to directly enter a group, 
restrict the possibility to leave to isolated agents, etc.). 
In this framework, two important points have to be con- 



• The decision for an agent to leave or enter the 
system can simply be given by constant probabil- 
ities. In the resulting dynamics, the number of 
agents who are present is simply fluctuating around 
a constant value that depends on these probabili- 
ties. Another possibility consists of fixing at each 
time step the number of present agents N(t). The 
imposed N(t) can be a given function of time such 
as for instance a periodic signal (possibly with su- 
perimposed stochastic noise) to mimic circadian 
rhythms. To mimic a realistic process, N(t) can 
also be given by an empirical time series from a 
real data set, as we now consider (note that the to- 
tal number TV of agents has to be at least equal to 
maxf TV(i)). 

In order to impose the number of agents present in 
the system at each time, we first construct the empirical 
time series. To this aim, various time steps can be used. 
We consider the natural temporal resolution obtained in 
SocioPatterns deployments, that is, 20s, but other reso- 
lutions could be chosen as well. Considering that agents 
act independently, and make choices in a simultaneous 
way, we then identify the empirical time step with one 
attempted status update per present agent. After each 
series of TV(t) attempted status updates (according to 
the rules described in Sec. IIII[) . we check if the num- 
ber of present agents is smaller or larger than the desired 
TV(t+l). If it is larger, we remove at random agents (i.e., 
put them in the "absent" state) in order to match the 
desired TV(t + 1). If it is smaller, absent agents are intro- 
duced into the system, and put into the isolated state. 
The system evolves in this way for a number of steps 
equal to the number of time steps of the empirical data 
set, and we monitor the time evolution of the number 
of contacts between agents, and of agents in each state, 
as a function of time. We also note that this model can 
in fact be simulated for arbitrarily long times by simply 
repeating the imposed time series N(t). 

Figure [5T] compares the resulting activity patterns with 
the empirically monitored activity, for two values of the 
parameter set (bo,bi,\), when TV(i) is taken from the 
data gathered at the scientific conference ESWC de- 
scribed in Sec. HU We note that, although only TV(t) 
is imposed to be exactly the same in the model and in 
the data, it is possible to tune the parameters so that 
the other measures (number of isolated nodes, of links, 
of triangles) are simultaneously similar to real data. Fig- 
ure illustrates that the system's dynamics is highly 
non-stationary. As shown in [22| . empirically observed 
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FIG. 21: (Color online) Timelines of the number of (from 
top to bottom): nodes and links in the instantaneous net- 
work (top), isolated nodes, groups of two nodes, groups of 
three nodes. The left column corresponds to the model with 
N(t) imposed from the ESWC dataset and two sets of param- 
eters, namely (bo,bi,A) — (0.55,0.8,0.9) (black curves) and 
(0.7, 0.7, 0.8) (red curves), the right column to the real ESWC 
data set. 



distributions of contact durations are however stationary. 
We check that this is also the case in our model by mea- 
suring the distributions of contact durations and of time 
spent by agents in each state in various time windows 
of different lengths, during which N(t) strongly varies. 
As shown in Figs. [22] and [23j for a given parameter set 
(bo,bi,X), these distributions are broad, as in the origi- 
nal model with constant number of agents, and do not 
depend strongly on the imposed N(t), and can be super- 
imposed from one time window to the next. As in the 
real data, the only differences come from different cutoffs 
stemming from different statistics in the different time 
windows. 

As for the study of empirical data (see Sec [TTJ) , we also 
construct and study the aggregated networks of contacts 
between agents, on different time windows. Figure I24| 
to be compared with Fig. [T5] of Sec. IIII1 illustrates the 
basic properties of these networks, which are very similar 
to the empirical ones shown in Fig. [3] of Sec. HI With 
respect to the case of constant N of Fig. [TH the degree 
distributions shift more slowly toward large degrees, and 
remain broader, so that the average strength of nodes of 
degree k, [s{k)\ and the average participation ratio of the 
nodes of degree k, [ky 2 (k)}, keep more realistic shapes. 

We finally illustrate the versatility of the modeling 
framework by considering the case in which an agent 
who leaves the network cannot re-enter it at future times. 
Such additional assumption can be considered to model 
environments in which there is a flux of individuals, such 
as a museum. As shown in [23[ , the daily aggregated net- 
work of the interactions between individuals is then not 
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FIG. 22: (Color online) Distributions of contact durations, 
for N(t) given by the ESWC data set, and for the model with 
constant N, for various time windows. Parameter values are 
(bo, bi, A) = (0.55,0.8,0.9). 
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FIG. 23: (Color online) Distributions of the time spent in 
various states, for N(t) given by the ESWC data set, and for 
the model with constant N. Parameter values are (bo, bi, A) = 
(0.55,0.8,0.9). 



a small-world, as visitors entering the museum at very 
different times do not encounter each other. Figure [25] 
exemplifies how this kind of configurations can also be re- 
produced by our modeling framework, by simply impos- 
ing the empirically measured number of visitors present 
on the premises, N(t), at each time step, and that an 
agent leaving the system does not re-enter it. The re- 
sulting network has an elongated shape whose topology 
is dictated by the timeline of the visits, exactly as in [231 ] . 
and is strongly different from the case in which agents can 
leave and re-enter the system, as in a conference. 
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FIG. 24: (Color online) Aggregated networks' characteristics 
for N(t) given by the ESWC time series, and 60 = bi = 0.7, 
A = 0.8, for various time window lengths. Top left: degree 
distributions; top right: links' weights distributions; Bottom 
left: average strength (s(k)) of nodes of degree k. Bottom 
right: average participation ratio {fcj/2 (fc)) of nodes of degree 
k. 



V. CONCLUSION 

In this paper, we have studied with analytical and nu- 
merical means an agent-based model aimed at describ- 
ing the dynamics of human interactions in social gather- 
ings, as can be measured by recent technological wear- 
able sensors that measure proximity patterns of indi- 
viduals. This model is based on simple mechanisms in 
order to be easily implcmentable and extended through 
the introduction of more realistic and involved ingredi- 
ents. The resulting distributions of contact durations and 
of inter-contact times can be either narrowly or broadly 
distributed. A detailed study of the latter case, obtained 
by a mechanism reminiscent of preferential attachment, 
reveals moreover interesting non-equilibrium transitions 
between stationary and non-stationary phases. 

We have illustrated the versatility of the model by in- 
troducing two variations. In the first one, agents have a 
priori different propensities to interact, in order to mimic 
the heterogeneity in the social behaviors of individuals. 
We have shown that the phenomenology of the model 
is then conserved, with broad distributions of contact 
times, inter-contact times, and lifetimes of groups. In 
the second extension of the model, the population size 
can fluctuate. In this case, the population size can be 
taken as input of the model, as given by an empirical time 
series coming from a real-world data set. We have then 
shown that the model produces non-stationary dynami- 
cal networks whose features are close to the empirically 
observed ones. 

The present modeling framework, and in particular its 
extension to a varying population, opens various perspec- 
tives. For instance, we have considered indistinguishable 




FIG. 25: (Color online) Example of aggregated network of 
157 nodes obtained by imposing the timeline of the number 
of agents present at each time during a given day of data 
gathered during a SocioPatterns deployment at the Science 
gallery in Dublin Here (bo, 61, A) = (0.55,0.8,0.9), and 

an agent who leaves the network cannot re-enter it. The nodes 
are colored and positioned according to the entry time of the 
corresponding visitor, as in [2^|. from red (top right) to green 
(right) to blue (left) to violet (top left). The elongated shape 
of the network is similar to the one observed empirically in 




agents, who enter and leave the system at random times. 
It would however be possible to impose for each agent the 
time intervals in which he/she is present, which could 
be taken from real data. Such a procedure would pro- 
duce artificial data sets that closely mimic the empirical 
timelines. Although this would happen at the cost of 
a rather large input of empirical information, it would 
produce data sets that retain the details of the presence 
properties of individual agents of the empirical data sets, 
but with tunable contact and inter-contact time distri- 
butions. Another interesting outcome of the model with 
varying population size is that it makes it possible, start- 
ing from an empirical data set that is by definition limited 
in time, to create a dynamical network on arbitrarily long 
timescales, with the same properties as the real one, by 
simply repeating the time-series N(i) as many times as 
required. This corresponds to an interesting way of creat- 
ing a non-stationary dynamical network, without having 
to repeat the real sequence of contacts. On each new rep- 
etition of N(t), the model will generate a new sequence 
of contacts. 

The modeling framework we have presented, and the 
above-mentioned perspectives, are particularly interest- 
ing in the perpective of the study of dynamical processes 
on dynamical networks. The ability to tune the networks' 
properties is indeed crucial to understanding how each of 
these properties affect the dynamical processes. Gener- 
ating artificial data sets that are based on empirical ones, 
preserve a certain number of their properties, modify oth- 
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crs in a tunable way, and can be extended to large sizes 
or long times, represents therefore a very important step 
in such studies. 
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Appendix A: Self-consistent solution of the pairwise 
model 



In this section we give the details of the self-consistent 
calculation that is able to solve for the mean-field dy- 
namics of the pairwise interaction model. As explained 
in the main text, the rate equations Eqs. ([5]) for this 
model are solved together with the definition of the tran- 
sition rates ttiq (t) and -kqi (t) given by Eqs. ([7]) by making 
the self-consistent assumption Eq.©. 

Inserting in the definition of ttiq (t) and ttqi (t) given by 
Eqs. © the structure of the solution of the mean- field 
dynamical Eq. and the self-consistent assumption 
Eqs.®, we get 



7T 10 (t) = STTOX-^I- 
t' = l v 



t-t' 

N 



-2bi-l 

(Al) 



where B is the (3 function. Inserting (|A4[) into (|A2[) we 
get 



For large N we can evaluate (|A1|) by going to the con- 
tinuous limit. Therefore in Eq. (|A1|) we substitute the 
sum over time steps t' with an integral over the variable 
y' = t' /N. The transition rate 7Tio(y) = iV7Tio(i), that 
is, the average number of agents that shift from state 
1 — > in the unit time y = t/N, can be evaluated by the 
following integral: 

7T 10 (y) = 2m Q1 b iy - a ^ 2bl f x~ a (l + y- 1 - x)- 2b ^ x dx 

Jo 



= 2Nn ib 1 y~ ai f(a 1 ,2b 1 + l : y), 
where /(a, b, y) is given by 



(A2) 



f(a,b,y)=y-^ x~ a (l + y' 1 - X y b dx. (A3) 



The asymptotic expansion of /(a, b, y) for y ^> 1 is given 
by 



/(a, b, y) = — — r + B(l — b,l — a)j/- b + o(- + y~ b 
b-l \y 



Kw(y) = NTT 01 y ai . 



(A5) 



This expression proves that the self consistent assump- 
tion given by Eq. ([8]) is valid. In particular since we have 
assumed 

7Tio(y) = Nn 10 y- ao 7roi{y) = A^oijT" 1 

these relations are consistent with the result of Eq. (|A5[) 
obtained in the limit N — > oo, y 2> 1 if 



(A6) 



In order to find the expression for a and tt we use the 
conservation of the total number of agents. Indeed we 
have 



E 



N (t,t') + N^t') 



= N 



(A7) 



Using the Eqs. ©, (|A5[) and (|A6[) and substituting in Eq. 
(|A7|) the sum over t' with an integral over the variable 
x = y' /y, we get, in the limit N ^> 1 



(2&o-i) / x -<*(i + y - x )- 



2b 



-y 



-(2bi-l) 



x~ a (l + y - x)- 2bl dx 



dx 



N ,(A8) 



which yields 

ny- a (f(a, 2b Q , y) + /(a, 2b u y)) = 1 (A9) 

Finally using the asymptotic expansion Eq. (|A4[) we get 
the solution given by the Eqs. © that we rewrite here 
for convenience 

a = max(0, 1 - 2&i, 1 - 2b ) , 
sin [2ir min (bO, 61) 



[l-*(«,0)] 



(2b - l)(26x - 1) 
2(b + 61-I) 



5(a,0). 



(A10) 



(A4) 



Appendix B: Self-consistent solution of the general 
model 



In this appendix we solve the general model in which 
groups of different size are allowed and the parameter 
A is arbitrary. The strategy that leads to the solution 
of the mean-field equation of this dynamics is essentially 
the same as in the pairwise model but a new phase tran- 
sition occurs when A < 0.5. The dynamical Eqs. (fTS"]) 
can be solved as a function of the variables TT mn (t) by 
Eqs. (TTJJ and Eqs. (fT()|) assuming self-consistently that 
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that e(t) = e in the large time limit. In order to find the 
analytic solution of the mean-field dynamics it therefore 
important to determine the relations between the tran- 
sition rates ir mn (t) and the variables N n (t,t'). These 
relations are given by 

m, (t) = 2A^ ^ V (M') 
v 

7r„ l0 (t) = "y V (M'), «>2 

t' 

7r n+ll „(t) - (n + l)A ^" + ^ M V (M'), »>1 



JV 



7T 



i(*) = 2^ 



JVo(M') 
N 



Po(t,0 



7ro,„(*) = (l-A)E^^ft(tO, ^>2 
t' 

7r„,„ +1 (i) = („+!)(!- A) X) 0, «>1- 



(Bl) 



The coupled Eqs. fT3]). pi)) and flBT) can be solved by 
making the additional self-consistent assumptions on the 
transition rates ir mn (t) given by 



7Tmn(2/) = N7T mn y a ™ n 



(B2) 



where y = t/N and ir mn {y) = Nir mn (t). 

Applying the same technique as in Appendix [3] we can 
prove that all the exponents a m , ra are equal and given by 
Oim,n = Performing straightforward calculations we 
get the following relations 

7T„,o(n + 1) = X[w n -i,n + TX n+ x <n + 7f ,„] for n > 2 

7^1,0 = A[7Ti. + 7T2,l] 

7T„,o = (1 - A)tt„_ 1i0 + A7r rl+ i i0 for n > 3 
1-A. 



-TTl.O + A7T 3 .o 



12,0 



7fi = A-7T1 n + 2A7T2,0- 



(B3) 



Therefore if the self-consistent assumption is valid, the 
number of agents N n (t,t') in state n since time t', is 
given at time t by 

N {t,t') = M i 



A" 



JV 



AW 



7ri,o(i')A , 



A 



JV 



-2bi 



A \ N 

where the variable A is defined by 



^,0 = ( " + ^° y) ri + ^v (w+1) V) 



(B5) 



A 



Using the relations given by Eqs. (|B3|) we find 

n-1 

/ 1 — A > 



tt,i,o = 7^1,0 ( — r~) for n > 2. (B6) 



A 

Substituting Eq. flB6j| in the definition of K, Eq. (|B5]l , 
we find that A is only defined for A > 0.5. For A < 0.5 the 
summation in Eq. (|B5[) is in fact divergent and there is 
a breakdown of the self-consistent assumption Eq. (|B2I) . 
For A > 0.5 we can perform the summation and we get 



A = 



2(2A-1) 

3A- 1 
1 



2A- 1 



(B7) 



Finally the value of a and 7x1,0 are found by enforcing the 
conservation law of the number of agent N 

t 

££ JV »(*> t, )=-M (B8) 

t> = l n 

Therefore, in the large y limit i/> 1 we get the solution 



3A — 1 

a = max ( 0, 1 — bo— -, 1 — 2bi 



'2A-1 

The value of 711,0 depends on the value assumed by a 
(1) For a = 0, the value of fr^o is given by 
1 



(B9) 



7Tl,0 



2(^0 - m) 



2A ^ 



n + 1 



(n+ l)6i - 1 V A 



1-A 



-1 
(B10) 



(2) For a = 1 — &o 1a— 1 ' tne vame °f ™l,o 1S given by 
2(2A-1) 1 



3A-1 B (l-b ^,b ^) 
where B(a, b) indicates the Beta function. 
(3) For a = 1 — 2&i, the value of tti.o is given by 

A 

Tin 



(Bll) 



(B12) 



5(1 - 2&i,26i) 
where -B(a, 6) indicates the Beta function. 
The average coordination number is defined by 

t TV 

< n > = EE nA ^')- ( B13 ) 

t' n=0 



Substituting Eqs. (|B4} to the definition of (n), Eq. (|B13|) 
and applying the same transformation in Eq. (|A2I) to 
evaluate the integral over t, we get 



n=l 



) = E x (n + 1 ~ ^y~ a ^ ( n + (B14) 



18 



where y = t/N and f(a, b, y) is defined in Eq. (|A3|) . Sub- 
stituting the asymptotic expansion Eq. (|A4[) into (|B14I) . 
we get 



N . 

TTnO 



1 



(n + l)6i 



+ B(l-(n + l)6 1 ,l-a)y 1 - (n+1)fcl 



where -B(a, 6) indicates the Beta function. In the asym- 
totic limit y -> oo, using Eqs. (jBo) . (|B10j)-f|B12| ) and 
counting only the leading terms in Eq. (|B15|) to com- 
pute (n) for different value of a, we can recover Eqs. 



Appendix C: Self-consistent solution of the 
heterogeneous model for A = 1 

In this appendix we show the self-consistent calcu- 
lations that solve analytically the heterogeneous model 
with pairwise interactions. 

We assume self-consistcntly that the transition rate 

Ttlo(t) and ttqi decay in time as a power-law, i.e. we 
assume 

<'(t) = A V A V '^' (±y aM) (CI) 

Inserting this self-consistent assumption and the struc- 
ture of the solution given by Eqs. (|25[) in Eqs. (j2l))) we 
can evaluate tt 11 ^ in the limit N — > oo. Therefore we 
get, 



=m -a( v ,v') 



2N 



v'tiy-^fiaWM + hy) (C2) 

where /(a, 6, y) is given by 

f(aAy)=y- [h - 1) I x- a {l + y- l -x)- b dx. (C3) 



The asymptotic expansion to f(a, b, y) for y 3> 1 is given 
by 

f(a, b, y) = -i- + B(l - b, 1 - a)^" 6 + O ( - + y~ b 
b-1 \y 

(C4) 



where B is the Beta function. Inserting (|C4[) into (|C2|) . 
we get in the limit y > 1 

N 



. -0(17,17') 



2C(y) 



^y-aoW^y-^V) (C5) 



Similarly, inserting (|C1[) into the definition of C(y) given 
by Eq. (|27l) we get, in the limit y > 1 



(C6) 



where we make use of the asymptotic expansion (|C4|) . 
In the limit t;> 1, the integral above can be calculated 
approximately by the saddle point method if 7r^ changes 
with 77 much slower than y" 01 ^. Therefore we have 



(B15) where 7 and 77* are given by 



7 = mino:^) 



argmin^a(r?). 



(C7) 



(C8) 



V 

By comparing both sides of Eq. (|C5[) and using Eq. (|C7 
we get 



"01 
ot(v,V r ) 



1 ~ 77 ~ T?' 



0(77) + a(r?') + 7. 



(C9) 



Finally, in order to fully solve the problem we impose the 
conservation laws of this heterogeneous model. In partic- 
ular the total number of agent with value rji G (77, 77 + A77) 
is given by the following relation, 

J2[No(t,t',r ] ) + Y^N 1 (t,t',r,,r 1 ')}=NA(r,). (CIO) 

V r,' 

Inserting the self-consistent anzatz Eq. (|C1[) for ttqi (t) 
and Eq. (|C5|) into Eq. (|C10|) we get, in the continuous 
limit approximation valid for JV»1, 



7r 10» 



0(277 — 1) 

2rj - 1 

xB(l-2?7,l-a(77))77 1 - 2 ''- 



2r?) 

-1 



(Cll) 



where 



N 



2C(y) 



(1 - 77 - 77') 



1 — 77 — 77' 



+ 0(77 + 7/ - 1) 



x B(77 + r/-l,l-a(77'))?/ r ' + "' _1 



x 4,2r«("'W- (C12) 

We compute I {if) defined in Eq. (IC12I) by counting the 
leading term only. Therefore we find 

a(rj) = max(0, 1 - 2rj, 77 - 1 + 7 + D) (C13) 

with D given by 

D = max[?7 - 0(77)]. (C14) 



Solving the Eqs. (|C13|) and (|C14I) we get 7 = and 
D = \ and 77* = 1/2. Therefore we can determine the 
exponent a(rf) and 01(77, 77') that are given by 



(C15) 



a(r/) = max ^1 — 2?7, 77 — 
0(77,77') = 0(77) +a(r]'). 
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Moreover the constants 7f^ are given, in the limit W>1 The only solution to the above expression is D 
and y 3> 1, by 



Similarly, for 7 + D > i, 



*10 = 



, eM n < - 

B(l-2ri,2rf) '1—2 



B(r,-i,l) 

Solving equation (|C18[) . let 7 
1 - 2rj 



D < i, then 



(C16) 



( 

'/ 



a(rj) 
and 

77 - a(r?) 

obviously, 7 
77 = 1 - 7 



i<r/<l-7-D (C17) 
1 + 7 + D 77 > 1 - 7 - -D 



3?? - 1 77 < i 

77 5<?7<l-7 
1 — 7 - D 77 > 1 - 7 - D 



D (C18) 



and D is reached either at 77 = i 
D, so 



or 



D = max(-,l- D) 



(C19) 



a{rf) = 



1-2t? V < 

77 — 1 + 7 + D ?7> 



, 2-7-D 

?^ 3— 

2-7-D 

3 



77 - a(7?) 



377-1 77 < 

1 - 7 - D n> l^zR 



Both 7 and D are reached at 77 = 2 ** D , so 



7 = !_ 2(2- 7 -£>) 
D = (2 — 7 — £)) — 1 



(C20) 



(C21) 



(C22) 
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